function res = ITCGguidanceloop(resp)
global hInt hGui N1 N2 ku kv kr kdelta Td T0
t = resp(1) + hInt;
statep = resp(2:7);
tPNGp = resp(10);
etp = resp(13);
x = statep(1);
y = statep(2);
z = statep(3);
V = statep(4);
fpa = statep(5);
azi = statep(6);
SMT = sqrt(x^2+y^2);
LMT = sqrt(x^2+y^2+z^2);

%% calculate LOS angle lambda and phi
lambda = atan(z/SMT);
if x >= 0
    phi = asin(y/SMT);
else
    phi = -(asin(y/SMT)+pi);
end
% look angle
delta = fpa + lambda; % 铅锤面速度前置角
eta = azi + phi + pi/2; % 偏航方向速度前置角

% if t > 35
%     kr = 1e-7;
% end
% if t > 30
%     kr = 1e-6;
% elseif t < 5
%     kr = 1e-5;
% else
%     kr = 5e-6;
% end
% kr = 1.25e-10*t^2 - 2.52e-7*t + 1e-5;
% calculate tPNG
statea = statep;
tPNG = tPNGp;
if isnan(etp)
    for j = 1:1e4
        stateb = PNGloop(statea);
        statea = stateb;
        if statea(3) <= 0
            tPNG = (j-1)*hInt;
            break;
        end
    end
else
    for j = 1:(tPNGp/hInt)
        stateb = PNGloop(statea);
        statea = stateb;
    end
    [r, dr] = calRandDR(statea);
    tPNG = tPNGp - (kr*r*dr + 1)*hGui;
end

% the tracking error et
et = tPNG - (TR(Td, T0, t) - t);
% et = tPNG - (Td - t);
% the gain adjustment term u 式17
N3 = N2 + (N2-1)*ku*et;

% lambdadot:俯仰视线角的变化率 phidot:偏航视线角的变化率
lambdadot = V*(sin(fpa)*cos(lambda)+cos(fpa)*sin(lambda)*cos(eta))/LMT;
phidot = V*cos(fpa)*cos(azi + phi)/LMT/cos(lambda);

% 纵向通道
% v = kv*(abs(et)-et)*exp(-kdelta*abs(delta))*sign(delta);
% fpadotd = -N3*lambdadot + v;
% azidotd = -N2*phidot;
% 侧向通道
v = kv*(abs(et)-et)*exp(-kdelta*abs(eta))*sign(eta);
fpadotd = -N1*lambdadot;
azidotd = -N3*phidot + v;

% calculate the guidance commands
[alphad, bankd] = calCommands(z, V, fpadotd, azidotd, fpa);
state = rk4(@TerminalDynamic, statep, alphad, bankd, hInt);
[r, dr] = calRandDR(state);

res(1) = t;
res(2:7) = state;
res(8) = lambda;
res(9) = phi;
res(10) = tPNG;
res(11) = alphad;
res(12) = bankd;
res(13) = et;
res(14) = r;
res(15) = dr;
res(16) = delta;
res(17) = eta;
end